Stratifications and foliations in phase portraits of gene network models

Periodic processes of gene network functioning are described with good precision by periodic trajectories (limit cycles) of multidimensional systems of kinetic-type differential equations. In the literature, such systems are often called dynamical, they are composed according to schemes of positive and negative feedback between components of these networks. The variables in these equations describe concentrations of these components as functions of time. In the preparation of numerical experiments with such mathematical models, it is useful to start with studies of qualitative behavior of ensembles of trajectories of the corresponding dynamical systems, in particular, to estimate the highest likelihood domain of the initial data, to solve inverse problems of parameter identification, to list the equilibrium points and their characteristics, to localize cycles in the phase portraits, to construct stratification of the phase portraits to subdomains with different qualities of trajectory behavior, etc. Such an à priori geometric analysis of the dynamical systems is quite analogous to the basic section “Investigation of functions and plot of their graphs” of Calculus, where the methods of qualitative studies of shapes of curves determined by equations are exposed. In the present paper, we construct ensembles of trajectories in phase portraits of some dynamical systems. These ensembles are 2-dimensional surfaces invariant with respect to shifts along the trajectories. This is analogous to classical construction in analytic mechanics, i. e. the level surfaces of motion integrals (energy, kinetic moment, etc.). Such surfaces compose foliations in phase portraits of dynamical systems of Hamiltonian mechanics. In contrast with this classical mechanical case, the foliations considered in this paper have singularities: all their leaves have a non-empty intersection, they contain limit cycles on their boundaries. Description of the phase portraits of these systems at the level of their stratifications, and that of ensembles of trajectories allows one to construct more realistic gene network models on the basis of methods of statistical physics and the theory of stochastic differential equations.


Introduction
At present time, investigation of questions of existence of periodic trajectories (cycles) in phase portraits of systems of non-linear differential equations simulating functioning of various natural processes is carried out in most fields of applied mathematics. Detection of such cycles, their localization in the phase portraits, description of their characteristics, such as stability, (non)uniqueness, etc. have a long history (Poincaré, 1892). These problems have generated a whole range of research directions in pure mathematics: qualitative theory of differential equations, theory of dynamics systems, etc., which in turn have a great impact on corresponding applied disciplines. At their junction, the famous 16-th Hilbert's problem, and the "center-focus" problem, related to seemingly just a pictorial case of two differential equations with two unknown functions of one variable (time) have appeared.
Here, in the present paper, we study systems of kinetic equations of higher dimensions, considered as functioning of circular gene networks models: It is assumed here and below that j = 1, 2, …, n; n ≥ 3, and that j -1 = n, if j = 1. In all these equations, non-negative functions x j (t) denote concentrations of species in the gene networks, and positive coefficients k j characterize the rates of their degradations (Likhoshvai et al., 2020).

Consider the system (1) in the vector form
where the vector-function X(t) is defined by its coordinate functions x j (t). The divergence of this vector-field F(X ) is constant and negative: div F(X ) ≡ -k 1 -k 2 -… -k n < 0.
It is well-known (Arnold, 1989) that in this case, n-dimensional volume of any finite domain in the phase portrait decreases exponentially during the shifts of its points along trajectories of the system (1) as t grows. This does not mean that each such domain collapses to a point. For the dynamical systems considered here, these limit sets are two-dimensional invariant surfaces in their n-dimensional phase portraits.
We call the dynamical system (1) block-linear if for all j each function f j which describes the rate of synthesis of the j-th component of the gene network is a step-function (threshold function) Here, a j are some positive constants. Decreasing functions L j describe negative feedbacks in the gene network and increasing functions Γ j correspond to positive feedbacks.
For one particular case k j = 1 for all j, investigation of cycles of similar block-linear systems was realized in (Glass, Pasternack, 1978;Akinshin et al., 2013;Ayupova, Golubyatnikov, 2014;Golubyatnikov, Gradov, 2021). Under the same assumptions, questions of existence of cycles in smooth analogues of these systems were studied in (Elowitz, Leibler, 2000;Glyzin et al., 2016;Kolesov et al., 2016) in the cases when these systems are symmetric with respect to cyclic permutations of pairs of the variables x j .
In recent publications Go lubyatnikov, Minushkina, 2019, 2020Likhoshvai et al., 2020;Ivanov, 2022), existence, uniqueness, and stability of the cycles of block-linear dynamical systems of some different dimensions with arbitrary positive coefficients k j were proved with the help of stratification of phase portraits to subdomains according to behavior of trajectories. It was shown there that these phase portraits contain cycles if and only if a j > 1 for all j and that the parallelepiped Q n = [0, a 1 ] × [0, a 2 ] × … × [0, a n ] in the positive octant of the space ℝ n is a positively invariant domain of the dynamical system (1). This means that trajectories of all points of the domain Q n do not leave it and that all cycles of the system (1) are contained in the interior of Q n . We consider below the dynamical systems of the type (1) in the case a j > 1 for all j only. Physical interpretation of this condition means that the maximal rate of synthesis of any component of the gene network exceeds that of its degradation.

Stratifications and foliations in phase portraits of gene network models
Here, each index ε j equals 0 or 1, and I j (0) = [0, 1], I j (1) = = (1, a j ]. Let E be the common point of all these blocks (all its coordinates equal one). In each of these blocks, the equations of the system (1) take the simplest linear form dx j dt = k j (x j -a j (1 -ε j-1 )), and solution to the Cauchy problem for this system has a simple representation In the present paper, for some low dimensional block-linear dynamical systems considered as models of gene networks functioning, we study the behavior of ensembles of their trajectories and show the existence of families of two-dimensional surfaces that are invariant with respect to shifts along trajectories of these systems and contain their cycles. This makes the qualitative analysis of trajectory behavior and interpretation of numerical experiments with these models much simpler.

Three-dimensional dynamical system
In the papers , we considered a 3D block-linear dynamic system: Trajectories of all points of the block {001} pass through six blocks of decomposition of the domain Q 3 from block to block according to arrows of the following diagram only: …→{001}→{011}→{010}→ {110}→{100}→{101}→{001}… (4) Denote by W 3 1 a union of blocks listed in the diagram, this is a positive invariant domain of the system (3), its interior is homeomorphic to torus. Note that trajectories of points of two blocks, {000} and {111}, eventually leave them in the invariant domain W 3 1 and further stay there. Thus, cycles of the system (3) do not intersect these two blocks . Stratification of phase portrait of the system (3) consists of two parts: the domain W 3 1 and the union of two blocks, {000}, {111}.
Consider a two-dimensional face F 0 = {001}∩{011} which separates the blocks {001} and {011} as well as other faces F m which separate incident blocks of the diagram (4): After transition along all six arrows of this diagram, trajectories of all points of the face F 0 return to it, each trajectory with its own time. Composition Ψ: F 0 →F 0 of all these six shifts from face F m to face F m+1 , m = 0, 1, 2, 3, 4, and F 5 →F 0 is called the Poincaré map.
The main technical result of the papers  is the following Lemma 1: а) the Poincaré map is monotonic: if for points A(v 1 ; v 2 ) and B(w 1 ; w 2 ) relations v 1 < w 1 and v 2 < w 2 , are satisfied then ψ 1 (v 1 ; v 2 ) < ψ 1 (w 1 ; w 2 ) and ψ 2 (v 1 ; v 2 ) < ψ 2 (w 1 ; w 2 ). For this partial order relation, we use a notation: A B, Ψ(A) Ψ(B); b) if w 1 and w 2 are sufficiently small then w 1 < ψ 1 (w 1 ; w 2 ) and w 2 < ψ 2 (w 1 ; w 2 ), i. e., B Ψ(B); c) at each point of the face F 0 , the first derivatives of the coordinate functions ψ 1 and ψ 2 are strictly positive and their second derivatives are strictly negative. This implies that the Poincaré map Ψ:F 0 →F 0 has two fixed points exactly; one of them is the point E 3 which lies at the boundary of F 0 and the other one, denoted by P * , is contained in the interior of the face F 0 . Trajectory of the point P * returns to this point after transition through the blocks of the diagram (4) and, therefore, it is a cycle. Since the map Ψ has just one nontrivial fixed point P * , the system (3) does not have any other cycles.
In the same paper, for the fixed points E 3 and P * of the Poincaré map, Jacobian matrices J 2 (E 3 ) and J 2 (P * ) were calculated and it was shown that the eigenvalues λ 1 (P * ), λ 2 (P * ) of the matrix J 2 (P * ) are different, positive and do not exceed one, which means exponential stability of the cycle of the system (3). We denote this cycle discovered in (Golubyatnikov, Ivanov, 2018) by ℂ 3 . Lemma 1 also implies that both these Jacobian matrices are positive, so it is possible to use the Frobenius-Perron theorem (Gantmacher, 1959) in our studies.
For sufficiently small ε > 0, we denote by T 2 ε ⸦ U(E 3 ) a triangle 0 ≤ u 1 + u 2 < ε with one vertex at the point E 3 and let F 0 be a truncated face F 0 \T 2 ε . Choose two segments [0, α 1 ] and [0, α 0 ] ⸦ [0, α 1 ] in this neighborhood so that α 1 = λ 1 (E 3 ) · α 0 . Let N 1 and N 0 , respectively, be the right endpoints of these segments, then Ψ([0, α 0 ]) = [0, α 1 ] and Ψ(N 0 ) = N 1 ; in the original coordinate system (w 1 ; w 2 ), the segments [0, α 0 ] and [0, α 1 ] are represented by arcs D 0 ⸦ D 1 with a common endpoint E 3 . Consider action of iterations of the Poincaré map to these arcs: The union D * of infinite sequence of mutually embedded arcs D k is a continuous monotonic arc connecting the points E 3 and P * ; after transition along arrows of the diagram (4), trajectories of points of D * return to this arc: the semi-interval D 1 \D 0 passes to semi-interval D 2 \D 1 which passes in turn to D 3 \D 2 , etc. Thus, trajectories of points of the arc D * generate an invariant (non-smooth) surface Σ 2 bounded by the cycle ℂ 3 761 СИСТЕМНАЯ КОМПЬЮТЕРНАЯ БИОЛОГИЯ / SYSTEMS COMPUTATIONAL BIOLOGY in the invariant domain W 3 1 ⸦ Q 3 . By the construction, this surface contains the point E 3 .
Starting such constructions of small segments [N 0 , N 1 ] in a neighborhood U(E 3 ) with points N 0 which do not lie on the axis E 3 u 1 and considering the images of these segments under iterations of the Poincaré map Ψ, we obtain a family of continuous monotonic arcs which leave the neighborhood U(E 3 ) and do not contain the point E 3 . For each pair of points N 0 , N 1 ⸦ U(E 3 )\ E 3 u 1 such that Ψ(N 0 ) = N 1 , the sequence N k = Ψ(N k-1 ) tends monotonically to the fixed point P * of the Poincaré map Ψ . Here, each segment [N 0 , N 1 ] generates, as above, a monotonic arc D * (N 0 ) being invariant with respect to the Poincaré map. Trajectories of points of such an arc, in their turn, form an invariant 2D surface Σ 2 (N 0 ) which intersects the surface Σ 2 by the cycle ℂ 3 exactly.
In a similar way, one can construct invariant surfaces which do not intersect the neighborhood U(E 3 ) in the domain W 3 1 . Let U(P * ) ⸦ F 0 be a neighborhood of the nontrivial fixed point P * , where the map Ψ can be linearized. We save the notations (u 1 ; u 2 ) for these linearized coordinates. For sufficiently small ε > 0, the Poincaré map transforms the ellipsis S 1 1 ⸦ U(P * ) with equation λ 1 (P * )u 2 1 + λ 2 (P * )u 2 2 = ε 2 to the circle S 1 0 with the equation u 2 1 + u 2 2 = ε 2 . Let I 1 (M 0 ) be a segment which joins the point M 1 S 1 1 with its image M 0 = Ψ(M 1 ) S 1 0 . All such segments are contained in U(P * ) in the ring between S 1 0 and S 1 1 . Each of these segments generates a sequence of continuous arcs D k (M 0 ), they are invariant with respect to the Poincaré map, and Ψ(D k (M 0 )) = D k-1 (M 0 ). For each of these arcs, trajectories of its points generate in W 3 1 an invariant surface bounded by the cycle ℂ 3 . Theorem 1. There exists two-dimensional invariant foliation in the invariant domain W 3 1 of the dynamical system (3); its leaves fill W 3 1 and contain the cycle ℂ 3 on their boundaries. One of these leaves contains the point E 3 .

Four-dimensional dynamical system
Recently, in the papers (Ayupova, Golubyatnikov, 2019;Golubyatnikov, Minushkina, 2021), we considered a fourdimensional block-linear system dx 1 dt In particular case, when k j = 1 for all j, questions of existence, uniqueness, and stability of cycles of such systems were studied in (Glass, Pasternack, 1978). Smooth analogues of similar systems were considered in (Hastings et al., 1977;Mallet-Paret, Smith, 1990).
An invariant domain Q 4 of the system (5) is decomposed by hyperplanes x j = 1 to 16 blocks {ε 1 ε 2 ε 3 ε 4 }. Blocks of this decomposition listed in the following diagram form an invariant subdomain W 4 1 in the phase portrait of (5) …→{1111}→{0111}→{0011}→{0001}→ {0000}→{1000}→{1100}→{1110}→{1111}→… (6) The arrows of this diagram show the only possible direction of trajectory transition from one block to another. The subdomain W 4 1 is one of two parts of stratification of the phase portrait of the system (5). For each block not listed here, trajectories of its points can pass to three adjacent blocks, two of them are contained in W 4 1 , and one is in Q 4 \W 4 1 (this is the second part of the stratification mentioned above). Algorithms of construction of such diagrams for the systems of arbitrary dimensions, both smooth and blocks-linear, are described in (Kazantsev, 2015;Kirillova, Minushkina, 2019).
As in previous sections, let us denote by 0 an intersection of two adjacent blocks {1111}∩{0111} in the diagram (6). After eight steps according to its arrows under shifts along trajectories, all points of this three-dimensional face return to 0 . Let Ψ 4 : 0 → 0 be a corresponding Poincaré map, T 3 ε ⸦ U(E 4 ) be a pyramid 0 ≤ u 1 + u 2 + u 3 < ε with the vertex at the point E 4 = (1; 1; 1; 1), and 0 be a truncated face 0 \ T 3 ε . In the paper (Golubyatnikov, Minushkina, 2021), it was shown that all statements of Lemma 1 are true for the map Ψ 4 , thus, this map has two fixed points exactly: E 4 and the point П * which is contained in the interior of the face 0 . This means that the invariant domain W 4 1 of the system (5) contains one cycle exactly, let us denote it by ℂ 4 . The following results were also established there.

Dynamical systems of higher dimensions
In the papers (Gaidov, Golubyatnikov, 2014;Ayupova, Golubyatnikov, 2021), we considered a five-dimensional block linear dynamical system 1 = L 1 (x 5 ) -k 1 x 1 ; 2 = L 2 (x 1 ) -k 2 x 2 ; … 5 = L 5 (x 4 ) -k 5 x 5 , (7) for which, as in previous sections, an invariant domain Q 5 = = [0, a 1 ] × [0, a 2 ] × … × [0, a 5 ] and its decomposition to blocks Stratifications and foliations in phase portraits of gene network models by the hyperplanes x j = 1 were constructed. Ten blocks of this decomposition form a stratum W 5 1 ⸦ Q 5 which is invariant with respect to shifts along trajectories of the system (7) passing through the blocks according to arrows of a cyclic diagram, similar to (4) and (6): {10110}→{10100}→{10101}→… Points of the four-dimensional face F 4 0 = {10101}∩{00101} under shifts along their trajectories after ten steps along the arrows of the diagram return to the face F 4 0 . For such a Poincaré map Ψ 5 :F 4 0 →F 4 0 , an analogue of Lemma 1 implies that the face F 4 0 contains two fixed points of this map exactly: the point E 5 = (1; 1; 1; 1; 1) and a point П 5 * in the interior of this face. The domain W 5 1 contains one cycle exactly. Let us denote it by ℂ 5 . This cycle is stable and passes through the point П 5 * (Ayupova, Golubyatnikov, 2021).
The magnitudes of eigenvalues of the matrix J 4 (П 5 * ) do not exceed one. In the case when these Jacobi matrices do not have any eigenvalues modulo equal to 1, construction of the invariant surface Σ 2 ⸦ W 5 1 with the boundary ℂ 5 and an invariant foliation in the domain W 5 1 is carried out exactly in the same way as above.
where b j are the maximum values of step functions Г j divided by the coefficients l j , j = 1, 2, 3, is decomposed by six hyperplanes m j = 1, p j = 1, j = 1, 2, 3, to 64 blocks which form a stratification of Q 6 to three subdomains, W 6 1 , W 6 3 , and W 6 5 , with different qualitative trajectory behavior.
The domain W 6 5 consists of 12 blocks, from which trajectories can transit to 5 adjacent blocks. In the symmetric case when k j = l j = 1, there are no cycles in this subdomain. However, W 6 5 contains a two-dimensional invariant surface consisting of piecewise linear trajectories attracting by the point E 6 = (1; 1; 1; 1; 1; 1) in a spiral way.
In the domain W 6 1 formed by 12 blocks, from which trajectories can enter one adjacent block only, the Poincaré map contains a unique non-trivial fixed point П 6 * , the trajectory of which is a stable limit cycle for all trajectories in this domain (Golubyatnikov, Minushkina, 2022).
In the domain W 6 3 which consists of 40 blocks, state transition diagram has a more complicated combinatorial structure. At present time, transitions of trajectories from one block to another in this subdomain have not been studied completely yet.
For smooth analogues of the dynamical system (8), the uniqueness of equilibrium point was established in (Ayupova et al., 2017). As in the case of block linear systems, hyperplanes passing through the equilibrium point and being parallel to coordinate ones decompose the invariant domain Q 6 to 64 blocks. If a linearization matrix of such smooth system in its equilibrium point has eigenvalues with positive and negative real parts and does not have any imaginary eigenvalues then the invariant domain W 6 1 contains a cycle ℂ 6 of this sytem (Ayupova et al., 2017). In the paper (Kirillova, 2020), conditions of existence of an invariant surface bounded by the cycle ℂ 6 in the domain W 6 1 were obtained.

Results of numerical experiments
The lefthand part of Figure shows 100 trajectories of the dynamical system (3). Each of these trajectories is contained in a corresponding leaf of the foliation in W 3 1 near the invariant surface Σ 2 . The values of parameters of this system are: k 1 = 0.4; k 2 = 0.3; k 3 = 0.6; a 1 = 1.3; a 2 = 1.4; a 3 = 1.7. The initial data are chosen in a random way in a rectangular neighborhood of the point It was shown in Ayupova, Golubyatnikov, 2021;Golubyatnikov, Minushkina, 2021;Minushkina, 2021) that trajectories of block-linear dynamical systems (3), (5), (7), (8) are piecewise smooth, the discontinuities of their derivatives are located on the planes x j = 1, this is clearly seen on the left part of Figure. In order to perform numerical simulations of trajectories of (3), we have developed a software project using the R programming language (https://www.r-project.org/) and the Shiny package (https://shiny.rstudio.com/). The source code is available on GitHub: https://github.com/AndreyAkinshin/pwLLL.
The simulations are performed in the cloud; the results are described at https://aakinshin.net/posts/dscs2/. The library ggplot (https://ggplot2.tidyverse.org/) is used here, as well as the package deSolve (http://desolve.r-forge.r-project.org/) that contains integration routines previously used to simulate other systems of gene networks. The user interface allows one to specify all parameters of the system (3).

Conclusion
In this paper, we have described a construction of invariant foliations, i. e. the families of invariant two-dimensional surfaces in phase portraits of low-dimensional block-linear models of circular gene networks. It was shown that on each leaf of СИСТЕМНАЯ КОМПЬЮТЕРНАЯ БИОЛОГИЯ / SYSTEMS COMPUTATIONAL BIOLOGY these foliations, trajectories of all its points are repelled by the boundary of the central part of the phase portrait and they are attracted by the limit cycle which describes an oscillating functioning of the corresponding gene network. Theorem 1 is illustrated by numerical experiments.
For the kinetic dynamical systems under consideration, the leaves of invariant foliations in the phase portraits play the role of level surfaces of collections of motion integrals studied in classical mechanics (Poincaré, 1892;Arnold, 1989). Reduction of dimensions of invariant subsets in the phase portraits allows us to give a digestible description of trajectories behavior and, in particular, simplifies considerably numerical experiments with such gene networks models (Likhoshvai et al., 2020). Construction of the foliations mentioned above and investigation of their geometric properties can be useful in studies of dynamical characteristics of more complicated models of gene networks functioning when a description of a big system is given on the basis of known results on its subsystems which have a simpler structure.